Strong surface contribution to the Nonlinear Meissner Effect 
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We demonstrate that in a d-wave superconductor the bulk nonlinear Meissner effect is dominated 
by a surface effect due to Andreev bound states at low temperatures. The contribution of this surface 
effect to the nonlinear response coefficient follows a 1/T 3 law with opposite sign compared to the 
bulk 1/T behavior. The cross-over from bulk dominated behavior to surface dominated behavior 
occurs at a temperature of T/T c ~ 1/y/K. We present an approximate analytical calculation, which 
supports our numerical calculations and provides a qualitative understanding of the effect. The 
effect can be probed by intermodulation distortion experiments. 

PACS numbers: 74.20.Rp, 74.25. N-, 74.45. +c 



In a superconductor with nodes in the gap function, 
quasi-particles near the gap nodes lead to an intrinsic 
nonlinear electromagnetic response In a d-wave su- 
perconductor this nonlinear Meissner effect appears as a 
linear magnetic field dependence of the magnetic penetra- 
tion depth at low temperatures 0, 0] , but can more sen- 
sitively be probed by temperature dependent intermodu- 
lation distortion or harmonic generation experiments Q . 
The nonlinear response coefficient shows an upturn at 
low temperatures following a 1/T law in a clean system 
down to temperatures of the order of 1/k, where k is the 
Ginzburg-Landau parameter of the superconductor. This 
behavior has been confirmed by intermodulation distor- 
tion experiments on high-T c cuprate superconductors 0- 
6]. At even lower temperatures nonlocal effects Q, as 
impurity effects Q , lead to a saturation of this low tem- 
perature upturn. 

So far, theoretical studies of the nonlinear Meissner ef- 
fect did not consider the special electronic structure that 
appears at the surface of a d-wave superconductor. At a 
surface that has a finite angle with the (100) direction of 
the crystal, Andreev bound states appear within a coher- 
ence length from the surface These states split in 
the presence of a screening current 13 - Tij and they carry 
an anomalous counter-flowing paramagnetic surface cur- 
rent 13, lij]. In previous work we have shown that the 
anomalous surface current leads to a strong modification 
of linear response properties 17 , 3] ■ Here, we study their 
influence on the nonlinear Meissner effect. We will show 
that the contribution of the surface Andreev bound states 
to the nonlinear response coefficient follows a 1/T 3 law, 
which will ultimately dominate the bulk 1 /T behavior at 
sufficiently low temperatures. We show that the cross- 
over from bulk dominated behavior to surface dominated 
behavior occurs at a comparatively high temperature of 
T/T c ~ 1/\/k. This means that even for a high n ~ 100 
as is realized in the cuprates the effect will become dom- 
inant at temperatures below about 0.1T C . 

In order to calculate the nonlinear response coefficient 
we solve Eilenberger's equations 3, 2(| fully momentum 



and energy dependent solving self-consistently the gap 
equation and the equation for the current density 



j(r) - AneN k B T V (v F (k)g(r, k, e„) ) (1) 

e„>0 FS 



together with the Maxwell equation 

V x V x A(r) = /xoj(r). 



(2) 



Here, A is the vector potential, Vf is the Fermi veloc- 
ity, Nq the single spin density of states, and g(r, k, e n ) 
the Eilenberger propagator on Matsubara frequencies e n . 
The full set of equations and a description of the numer- 
ical solution procedure based on the Riccati technique 
2l[ can be found in Ref. [H 

We consider a homogeneous superconducting half- 
space in the region x > with an external magnetic 
field Bo parallel to the z-axis, which shall be aligned 
with the c-axis of the crystal structure. In this geometry 
the current flows along the y-direction in the supercon- 
ductor. The gap function is assumed to have a rotated 
d-wave form A(x, 9) = Ao(x) cos2 (9 — a), where a is 
the angle of rotation with respect to the surface and the 
angle 9 denotes the direction of momentum within the ab- 
plane. As this problem is translationally invariant in y- 
and z-direction, all quantities only depend on the spatial 
variable x. The self-consistent solution of Eilenberger's 
equations on real frequencies allows us to calculate the 
local, angular resolved normalized density of states 



N (E, x, 9) = -Im g (x, 9, ie n -> E + i0 H 



(3) 



The equation for the y-component of the current den- 
sity |T]) can be transformed by contour integration and 
analytic continuation to the real axis 0, Q : 



dE / dO sin* 



j (x) = -eN v F 



f(E)[N + (E,x,9)-N_(E,x,9)} (4) 

where / (E) = 1+e 1 E / T is the Fermi function and N± de- 
notes the normalized density of states for comoving and 
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FIG. 1: Spatial dependence of the nonlinear coefficient 773 at 
T — 0.1T C and a = 7r/4 as a function of the distance a: from 
the surface in units of the coherence length £o- Inset: larger 
scale for x > 4£o, highlighting the sign change of 773. 



countermoving quasiparticles relative to the condensate 
flow, i.e. N + (E, x, 9) = A_ (E, x,9-ir). Once the cur- 
rent density distribution j(x) is obtained, the magnetic 
field distribution B{x) and the vector potential A(x) are 
found from integration of Eq. @. 

For a high-K superconductor the length scale of vari- 
ation of the vector potential, the magnetic penetration 
length A, is a factor of k larger than the variation of the 
Eilenberger propagator g on the length scale of the co- 
herence length £0 — hyp/ir Aq. Thus, for temperatures 
T/T c > 1/k it is a very good approximation to evalu- 
ate the angular resolved local density of states by a local 
Doppler shift of the fully nonlocal Eilenberger propagator 
in the absence of a vector potential, i.e. 

N±(E,x,0) = N(E±ev F - A(x),x,9) (5) 

where N (E, x, 9) on the right hand side is calculated with 
A(x) = but fully includes the surface Andreev bound 
states. Here, we have chosen the real gauge in which the 
vector potential is directly proportional to the superfluid 
velocity v s (x) = -^A(x). 

In order to determine the lowest order nonlinear re- 
sponse, Eq. ([5]) is substituted into Eq. (j4j and we make a 
Taylor series expansion of j in the vector potential A(x): 

-2e 2 v 2 F N a 771 (x) A(x) + 



(6) 



2e A v%N 

A 2 
^0 



773 (x) A 3 (x)+0(A 5 ) 



Here, the even terms in A cancel out due to symmetry. 
After a partial integration the dimensionless expansion 
coefficients are given by the expressions: 



1 + 



„ A 2 - 

0, 



dO sin 2 9 



d9 sin 4 9 



dE^N(E,x,9) (7) 



OO 

dE^N(E,x,9) (8) 



where Ao is the zero temperature gap value in the bulk. 
Note, that in contrast to the bulk calculation 0, Q the 
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FIG. 2: (Color online) Double logarithmic plot of 1 973 1 as a 
function of temperature T/T c for three selected positions: (a) 
x = 0, (b) x — 8£o, and (c) x = 45£o- The red dashed lines 
show a 1/T 3 dependence and the blue dotted lines a 1/T 
dependence. 



expansion coefficients now depend on the distance from 
the surface. Within a distance of the order of the co- 
herence length they contain contributions from the An- 
dreev bound states. The coefficient 771 describes the lin- 
ear response and the coefficient 773 the lowest order non- 
linear response. The spatial dependence of 773 is shown in 
Fig. [T]for a (110) surface (a = 7r/4) at a temperature of 
T = 0.1T C . Deep in the bulk, 773 is positive and reaches 
the low temperature value Aq/2T known from previous 
work However, when the surface is approached within 
a few coherence lengths, 773 changes sign and reaches ex- 
tremely large negative values at the surface. 

The temperature dependence of the modulus 1 773 1 is 
shown in Fig. Hon a double logarithmic scale for three 
selected spacial positions. Fig. [2jc) shows the temper- 
ature dependence at x = 45£o in the bulk. As is well 
known from previous work, 1 773 1 follows a 1/T law at low 
temperatures (blue dotted line). Right at the surface 
(x — 0), however, Fig. [^a) demonstrates that 1 773 1 is fol- 
lowing a 1/T 3 behavior (red dashed line). In Fig.^b) an 
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intermediate position at x = 8£o is shown. In this 
higher temperatures a 1/T law is followed. At a certain 
temperature, 773 changes sign and starts to follow a 1/T 3 
behavior below that temperature. These results clearly 
show that the nonlinear response coming from the sur- 
face area, where the Andreev bound states are present, 
is much stronger and of opposite sign than the nonlinear 
response in the bulk. 

In a typical intermodulation experiment only the total 
response of the system is probed. The quantity that is 
observed is the nonlinear change of the total inductance 
of the system The total inductance L can be calcu- 
lated from the total kinetic and magnetic field energy in 
the system via the equation 



-LI 



1 

2^o 



dx(B 2 (x)-n j(x)A(x)) (9) 



where / = dx j (x) is the total current per unit length 
0. Using Eq. ©, B = dA/d x, and the fact that 
the magnetic field vanishes in the bulk, Eq. (j9]) can be 
brought by partial integration into the more convenient 
form 



4° 
1 



(10) 



with Aq — A(x — 0). To lowest order in Aq the total 
current I generally will be of the form 



I = aiA + azAl 



(11) 



The intermodulation response is proportional to the non- 
linear coefficient 1^ A straightforward calcula- 
tion shows that this quantity can be related to the ex- 
pansion coefficients a\ and 03 using Eq. (|10|) 



d 2 L 
dl 2 



2a 3 



(12) 



1=0 



We have determined a\ and 03 from our numerical so- 
lution of Eilenberger's equations. The resulting values 

for are shown in Fig. H as a function of re- 

duced temperature for k = 63 (solid black circles) and 
k = 1000 (solid red squares) on a double logarithmic 
scale. Decreasing the temperature from T c , for a = 7r/4 
the nonlinear coefficient initially decreases and changes 
sign at a temperature near T/T c sa 2.4/^/k. Below that 
temperature the nonlinear coefficient increases following 
a 1/T 3 law and finally diverges at a temperature near 
T/T c rj 1/k. For comparison also the behavior for k = 63 
and a = is shown, when the surface states are absent 
(open circles). In this case there is no sign change and 
the nonlinear coefficient follows a 1/T behavior at low 
temperatures, as known from the bulk. 

In order to check the validity of the numerical calcula- 
tions and obtain a physical understanding of the results, 
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FIG. 3: (Color online) Double logarithmic plot of 

I UL 1 1=0 

as a function of temperature T/T c for « = 63 and a = n/4 
(solid black circles), for k = 1000 and a = 7r/4 (solid red 
squares), and for k = 63 and a — (open blue circles). The 
dashed lines show a 1/T 3 behavior and the dotted line a 1/T 
behavior. The solid lines show the approximation Eq. (|16p 
for k = 63 (black) and k = 1000 (red) , respectively. 



we made an approximate analytical solution of the prob- 
lem which we present now. For a piecewise constant gap 
function Eilenberger's equations can be solved analyti- 
cally [2ll ]. As an approximation we assume that the d- 
wave gap is constant in space. Then, the analytical solu- 
tion of Eilenberger's equations allows us to determine the 
residue of the zero energy pole of the Eilenberger propa- 
gator analytically, which contains the contributions from 
the zero energy bound states at the surface. As a result, 
we find the following expression for the bound state con- 
tribution to the local, angular resolved density of states 
for a — 7r/4 in the absence of an external field: 



N bs (E,x,0) =7rA |sin2<9|e~ 



^S(E) (13) 



The 5-function shows that the bound states are only 
present at zero energy. The exponential factor drops off 
on a length scale of the coherence length, showing that 
these states are localized at the surface. Introducing this 
expression into Eq. (|8]) the energy integration immedi- 
ately shows that the bound states lead to a 1 /T 3 scaling, 
which is of opposite sign than the bulk behavior, because 
^4 is positive at zero energy, but negative at higher 
energies. 

In order to determine the coefficients a\ and 0,3 in 
Eq. (fTTj) . we integrate Eq. (jSJ) using the following ap- 
proximations. For k = A/£o 3> 1 we can assume that the 
vector potential varies exponentially on the length scale 
of the penetration length A, and make the ansatz 



A(x) = (A - e)e- x / x + e e - 3x/x 



(14) 



The functions r\\ and 773 both vary on the length scale 
of the coherence length, which is much smaller than A. 
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Therefore, we can approximate them as 



c\5(x) 
036 (x) 



r]3b 



Here, 771b and 7735 are the bulk values of 771 and 773, re- 
spectively. The coefficients C\ and C3 describe the contri- 
butions of the surface bound states. They are obtained 
by substituting Eq. (TJ3]) into Eq. © and Eq. © and in- 
tegrating over x from to 00. This yields c\ — ~^fr-£o 



and C3 = — £q . The parameter e in Eq. (fT4"|) is dc- 



20T 3 



termined from the differential equation Eq. ^ together 
with Eq. © and up to order Aq found to be 



8 A§77i 6 



(15) 



With these approximations we find from integrated 
Eq. @ 

a\ = -2e 2 VpN {c\ + \r] lb ) 



^0 



A 

C3 + jr?3& 



Using the low temperature limiting expressions rjn, ~ 1 
and 773b ~ A /2T finally leads to 



d 2 L 

a/ 2 
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J=0 



16e 4 7j4 TV^AgA 3 



1 _ E. A An 

1 a .. rp 



(16) 



This expression shows that upon lowering the tempera- 
ture from T c the total nonlinear response of the system 
initially follows the 1/T increase caused by the bulk non- 
linearities (first term in the numerator). At a tempera- 
ture of the order of T/T c ~ 1/v^ the nonlinearities of 
the surface states become comparable with the bulk con- 
tributions and cancel them (second term in the numer- 
ator). Below that temperature the 1/T 3 increase with 
opposite sign dominates due to the surface states. Fi- 
nally, at a temperature of the order of T/T c ~ 1/k the 
nonlinear response diverges (denominator). This diver- 
gence signals the breakdown of the large k approximation 
we have used here. The approximate expression Eq. (|16[) 
is shown in Fig. [3] together with the numerical results. 
The agreement is quite good at low temperatures despite 
the approximations made. 

To conclude, we have shown that in a <i-wave supercon- 
ductor surface Andreev bound states lead to a strong con- 
tribution to the nonlinear Meissner effect, which follows 
a 1 /T 3 behavior at low temperatures and is of opposite 
sign compared to the bulk nonlinear response. At tem- 
peratures below T/T c r~j 1/ ' yfH these contributions dom- 
inate the total nonlinear response. Such temperatures 



are readily available in intermodulation experiments and 
make them a tool to study surface Andreev bound states. 
The fingerprint of the Andreev bound states should be a 
1/T 3 temperature dependence and a sign change (180° 
relative phase change) in the nonlinear part of the induc- 
tance. So far, intermodulation experiments have been 
mostly done on systems with (100) surfaces, where An- 
dreev bound states are absent. In systems with (110) sur- 
faces the effect studied here should become most promi- 
nent. 
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